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1.  Introduction 


Titanium  alloyed  with  6%  aluminum  and  4%  vanadium  (Ti-6A1-4V)  is  a  material  of  increasing 
importance  in  armor  applications  for  the  U.S.  Army  because  of  its  light  weight  (compared  to  steel) 
and  its  high  strength.  The  U.S.  Army  Research  Laboratory  (ARL)  has  developed  a  low  cost 
alternative  to  the  common  (but  expensive)  aircraft-grade  Ti-6A1-4V.  Extensive  Hopkinson  Bar 
tests  have  been  conducted  on  this  low  cost  material  (subsequently  referred  to  as  simply  Ti-6A1- 
4V),  including  Casern’s  recent  experiments  at  strain  rates  as  great  as  50,000/s  (7),  in  an  effort  to 
understand  the  material’s  dynamic  behavior.  The  ultimate  goal  is  to  develop  a  capability  to  predict 
the  material’s  strength  response  during  the  extreme  loading  conditions  encountered  during  ballistic 
impact  and  penetration.  This  goal  would  be  achieved  if  a  suitable  strength  model  were  established 
for  use  in  finite  element  and  finite  difference  shock  physics  codes,  such  as  EPIC1  (2)  and  CTH2 
(5). 

Toward  this  end,  several  existing  strength  models  were  considered  as  candidates  for  calibration 
with  the  use  of  the  new  Hopkinson  Bar  data:  the  Johnson-Cook  strength  model  and  the  Zerilli- 
Armstrong  series  of  models.  The  Johnson-Cook  model  ( 4 )  calculates  yield  strength  by  the 
relation 

T  =  [j  +  5f"][l  +  Cln/][l-r*m]  (1) 

In  equation  1 ,  A,  B,  C,  n,  and  m  are  material  constants,  e  is  the  strain,  £  *  is  the  non-dimensional 
strain  rate,  and  T*  is  the  homologous  temperature.  This  model  was  found  to  be  inadequate  for 
the  Ti-6A1-4V  material  since  it  does  not  model  the  nonlinear  strength  behavior  (versus  the 
logarithm  of  strain  rate)  exhibited  by  the  material  and  does  a  poor  job  of  modeling  two  material 
characteristics  (discussed  in  more  detail  later)  that  are  important  to  predicting  localization:  the 
coupling  of  thermal  and  rate  effects  and  stress  saturation  (approaching  an  asymptotic  value  or 
reaching  a  maximum  value)  at  high  strains3.  Nevertheless,  this  model  was  calibrated  with 
Casern’s  data  and  tested  in  simulations  of  penetration  experiments  as  a  reference  for  comparison. 
Details  are  discussed  in  a  later  section. 

The  Zerilli-Armstrong  models  are  physically  based,  and  there  are  several  generations  of  models. 
Initially,  the  model  addressed  metals  with  either  fee  (face-centered  cubic)  or  bcc  (body-centered 
cubic)  crystal  structures  (5,6).  The  yield  strength  is 

Y  =  A +  [ci+C2yfs]e{~Cl+c^}T  +C5sn  (2) 


^PIC,  which  stands  for  Elastic-Plastic  Impact  Code,  uses  a  Lagrangian  frame  to  solve  the  equations  of  motion. 
2CTH,  which  is  not  an  acronym,  uses  an  Eulerian  frame  to  solve  the  equations  of  motion.  Both  EPIC  and  CTH 
are  widely  used  across  the  defense  research  and  development  community  to  model  problems  in  shock  wave 
propagation. 

3Unless  otherwise  noted,  the  term  strain  refers  to  equivalent  plastic  strain. 
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in  which  s  is  strain,  s  is  strain  rate,  and  T  is  temperature.  By  appropriate  choice  of  the  constants 
(A,  Ci,  C2,  C3,  C4,  and  C5),  the  model  is  applied  to  either  fee  (C\  =C5=0)  or  bee  (C2=0)  metals. 

This  model  was  used  with  success  in  earlier  work  (7).  Zerilli  and  Armstrong  expanded  the 
applicability  of  the  fcc/bcc  model  by  the  development  of  a  newer  representation  for  hep  (hexagonal 
closely  packed)  metals  (8,9).  Then  they  extended  their  hep  model  to  address  shear  instability  (10), 
an  important  consideration  for  Ti-6A1-4V.  This  more  recent  hep  form  was  desired  here  for  two 
reasons.  First,  the  hep  alpha  phase  in  Ti-6A1-4V  dominates  the  microstructure  until  transformation 
to  the  bcc  beta  phase  occurs  at  approximately  996  °C  (11).  Thus,  hep  behavior  is  an  important 
component  of  the  material’s  ballistic  behavior,  and  hep  behavior  can  be  only  approximately 
described  by  the  fcc/bcc  model.  Second,  a  model  that  would  address  shear  instability  and  the 
resultant  strain  localization  was  desired  for  the  present  effort.  This  model,  calibrated  with  the  new, 
higher  strain  rate  data,  was  expected  to  give  better  results  than  were  obtained  in  the  earlier  work. 


2.  The  Modified  Zerilli-Armstrong  Model 


The  localizing  strength  model  proposed  by  Zerilli  and  Armstrong  (10)  for  the  strength  of 
predominantly  hep  Ti-6A1-4V  is  somewhat  different  from  equation  2: 

Y  =  C0  +  C,epT  +  C2  ^sr{l-e-£,£')  (3) 

Here,  Co,  Ci,  and  C2  are  material  constants,  a  and  p  are  functions  of  strain  rate  (described  next), 
and  T  is  the  absolute  temperature.  The  expression  under  the  radical  will  be  referred  to  as  the 
“strain  function”.  The  strain  is  s,  and  the  material  constant  sr  (named  the  “recovery  strain”  by 
Zerilli  and  Armstrong),  affects  the  strain  at  which  saturation  of  the  stress  is  achieved. 

The  square  root  comes  from  theoretical  (i.e.,  ideal)  considerations  of  dislocation  motion  during 
the  plastic  deformation  of  crystals  given  by  Taylor  (12).  In  that  work,  Taylor  concluded  that  the 
yield  stress  depends  on  the  square  root  of  strain  (and  other  factors).  The  problem  with  a  simple 
square  root  function  is  that  the  stress  will  not  saturate  at  high  strains  in  a  function  of  the  form 
Y  =  K\fs  .  This  type  of  strain  hardening  behavior  is  detrimental  if  modeling  of  localization  is 
desired,  since  localization  depends  on  thermal  softening  effects  overtaking  the  strain  hardening 
effects.  In  fact,  Ti-6A1-4V  does  saturate  at  high  strains,  so  Zerilli  and  Armstrong  replaced  the 
simple  strain  term  in  the  Taylor  strain  hardening  of  earlier  models  with  the  strain  function  in 
equation  3,  which  shows  saturation  at  high  strains. 

The  thermal  and  rate  effects  are  coupled  in  the  exponential  terms  in  equation  3;  the  parameters  a 
and  P  are 

a  =  a0  -ax  In  s 


P  =  Po~PMe 
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(4) 

(5) 


where  ao,  ai,  J3o,  Pi  are  material  constants.  Zerilli  and  Armstrong  state  that  s  is  the  strain  rate,  a 
dimensional  quantity.  Taking  the  logarithm  of  a  dimensional  quantity  is  at  best  awkward,  for 
example,  requiring  awkward  units  for  the  constant  pi.  Johnson  and  Cook’s  approach  to  the 
logarithmic  function  was  to  non-dimensionalize  the  strain  rate  by  a  reference  strain  rate  (1/s). 
Several  unsuccessful  attempts  were  made  here  to  include  the  same  approach,  but  none  of  the 
attempts  were  able  to  accurately  model  the  data  across  the  wide  spectrum  of  strain  rates 
(although  each  attempt  was  accurate  near  the  reference  strain  rate  chosen).  Finally,  the  following 
approach  was  adopted  (the  term  sp  is  borrowed  from  Zerilli  and  Armstrong  in  reference  8). 

This  approach  does  not  functionally  change  equations  4  and  5  but  only  substitutes  a  new  material 
constant  for  one  of  the  original  ones.  Let 

(«) 

In  sp 


where  ep  is  a  material  constant,  but  since  it  replaces  pi  in  equation  5,  the  model  (equation  3) 

retains  the  same  number  of  material  constants.  The  units  of  Pi  in  equation  6  are  unchanged  from 
the  original  form.  Substitute  equation  6  into  equation  5: 


0  =  0o 


In  s 

In 


(7) 


Similar  reasoning  yields  the  new  form  for  a: 


a  =  a0 


In  s 

ln4 


(8) 


The  awkward  units  of  the  material  constants  have  been  eliminated4.  For  present  purposes,  strain 
rate  units  will  be  s'1.  In  practice,  this  new  form  was  better  able  than  the  original  form  to  capture 
the  strain  rate  effect  on  the  strain  hardening  for  Casern’s  data  set.  Importantly,  the  constants  ea 

and  kp  have  physical  significance.  When  s  is  less  than  sa ,  a  is  positive  and  the  exponential 

term  in  equation  3  decreases  with  increasing  temperature.  This  “thermal  softening”  (reduction  of 
yield  strength  with  increasing  temperature)  leads  to  increased  strain,  which  leads  to  higher 
temperatures.  This  cascade  effect  may  lead  to  localization.  However,  if  s  is  greater  than  sa ,  a 

is  negative  and  the  exponential  term  in  equation  3  increases  with  increasing  temperature,  and  the 
cascade  does  not  occur.  A  similar  argument  applies  to  sp  .  Thus,  large  values  of  the  two 

(relative  to  s )  are  desired  for  modeling  strain  localization. 

For  Casern’s  data,  the  square  root  in  equation  3  provided  only  a  fair  fit.  Since  the  empirical  fits 
for  the  Johnson-Cook  model,  which  are  found  in  the  literature  for  various  materials,  give  a 


4Solving  equation  6  for  p0  and  substituting  that  into  equation  5  will  eliminate  the  logarithm  of  a  dimensional 
quantity,  but  the  fits  obtained  with  such  a  form  were  not  as  good  as  those  obtained  by  the  method  just  described. 
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variety  of  values  for  the  exponent  of  the  strain  hardening  term  (n  in  equation  1),  optimization 
was  conducted  with  the  radical  replaced  by  a  variable  exponent  (i.e.,  an  additional  material 
constant).  Surprisingly,  the  constant  was  repeatedly  precisely  1  in  a  variety  of  optimization 
runs5.  The  optimization  technique  is  discussed  in  the  next  section.  Analysis  of  equation  3  with 
an  exponent  of  1  shows  that  the  algebraic  form  still  saturates  at  high  strain.  In  fact,  the  model’s 
ability  to  capture  the  saturation  of  the  data  at  high  strain  rates  was  improved  over  the  fit  obtained 
with  the  square  root  function.  This  change  also  improved  the  representation  of  the  strain  rate 
effect.  The  overall  fit  to  the  data  was  greatly  improved,  and  this  is  quantified  later.  Therefore, 
the  square  root  in  equation  3  is  deleted  for  the  present  work. 

Thus  for  Casern’s  data  set,  the  modified  form  of  the  Zerilli- Armstrong  hep  strength  model  is 

Y  =  C0  +  CyeffT  +  C2  \sr  (l  -  es,e' )]  (9) 

with  p  and  a  described  by  equations  7  and  8.  This  function  fit  Casern’s  data  quite  well  and  is  the 
subject  of  the  remainder  of  this  report. 


3.  Calibration  of  the  Model 


The  initial  (plastic  strain  =  0)  yield  strength  of  the  material  at  T  =  OK  is  Yc.  Under  these  conditions, 
equation  9  reduces  to 

r=c0+c1  (10) 

Additional  unpublished  data  describing  initial  yield  behavior  of  the  Ti-6A1-4V  at  various  tempera¬ 
tures  from  8 1  K  to  693  K  were  extrapolated  to  0  K  with  the  use  of  an  exponential  function  of  the 
form  of  the  second  term  in  equation  9.  The  value  of  Yc  so  obtained  was  1356  MPa.  Thus,  a 
constraint  on  the  fit  of  equation  9  is  that 

C0+C1  =  1356  MPa  (11) 

Co  itself  has  physical  significance;  it  is  the  static  yield  strength  of  the  material  (i.e.,  s  — »  0  in 
equation  9).  However,  in  the  present  work,  Co  was  not  obtained  by  quasi-static  yield  tests  but 
was  fit  to  the  data  in  the  global  process  described  next. 

Casern’s  data,  summarized  in  the  next  section,  consisted  of  stress-strain  curves  for  five  different 
strain  rates6  from  960/s  to  46,400/s.  Each  of  the  five  data  sets  was  itself  an  average  of  several 
Hopkinson  Bar  tests.  For  present  purposes,  each  data  set  was  represented  by  a  discretization  of 


5For  example,  omitting  one  of  the  data  sets  (e.g.,  9585/s)  or  fitting  the  initial  yield  parameters  independently  before 
fitting  the  strain  hardening  parameters,  or  adding  other  material  constants,  or  trying  different  forms  for  a  and  p. 

6A  sixth  data  set  (0.1/s)  was  available  but  was  not  used  in  the  calibration  process  because  this  rate  is  not  of  primary 
importance  in  ballistic  problems,  and  including  it  would  have  worsened  the  fit  at  the  important  high  strain  rates. 
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the  data  set  into  roughly  20  stress-strain  pairs  taken  along  the  curve  between  initial  yield  and 
failure.  An  optimum  fit  of  equation  9  to  this  composite  data  set  was  sought.  The  Global  Local 
Optimizer  (GLO)  (73)  was  used  for  this  purpose.  Through  sophisticated  mathematics,  GLO 
varies  all  the  parameters  of  equation  9  (with  the  exception  of  C\,  which  was  eliminated  from  the 
optimization  process  through  use  of  equation  1 1)  in  a  systematic  way  and  optimizes  the 
predictions  of  the  resulting  equation  relative  to  the  composite  data  set  by  searching  for  a 
minimum  of  a  user-defined  figure  of  merit  (FOM).  Here,  the  FOM  was  defined  as 


FOM 


— y — y 

AT  ^  ^ 


® expt  ® calc 


N  M  cr 

iy  i  1V±  J  ^  expt 


(12) 


in  which  aexpt  is  the  experimental  stress,  acaic  is  the  stress  calculated  from  equation  9,  the  inner 
summation  (j)  is  over  the  20  stress-strain  pairs  in  each  strain  rate  data  set,  and  the  outer 
summation  (i)  is  over  the  N  =  5  strain  rate  data  sets.  This  FOM  (x  100)  is  actually  the  average 
percentage  error  between  the  experimental  stress  and  the  calculated  stress  at  a  strain  point.  The 
FOM  is  an  indicator  of  the  quality  of  the  fit  to  the  Hopkinson  Bar  data  and  is  not  necessarily 
indicative  of  the  model’s  performance  in  simulations  of  ballistic  events.  The  FOM  is  an 
indicator  of  the  quality  of  fit  when  all  strains  and  strain  rates  are  considered,  whereas  the 
important  part  of  a  ballistic  event  may  include  only  a  subset  of  those  (that  subset  taken  alone 
may  have  a  different  FOM).  A  larger  consideration  is  the  degree  to  which  Hopkinson  Bar 
experiments  model  ballistic  events. 


The  GLO  fitting  process  consists  of  selecting  a  set  of  material  constants  and  then  marching 
through  each  data  set  from  initial  yield  to  failure,  calculating  the  stress  at  each  experimental 
strain.  The  temperature  increase  in  moving  from  one  strain  to  the  next  is  required  during  this 
process;  it  is  assumed  to  be  entirely  attributable  to  adiabatic  heating  from  plastic  work  and  is 
calculated  from 

AT—2-cr, (13) 
pc„ 


where  i  denotes  the  current  stress-strain  pair  and  i-1  the  previous  pair.  This  implicit  scheme  is 
necessary  because  the  temperature  must  be  calculated  before  the  current  stress  (a;)  can  be 
calculated.  A  second  order  scheme  (estimating  a,  from  a,./  and  G;_2 )  was  found  to  be  unnecessary. 
For  the  present  Ti-6A1-4V  work,  coefficient  p  =  0.9,  p  =  4424  kg/m3,  and  cp  =  590  J/Kg.K. 
Adiabatic  heating  was  included  for  all  dynamic  strain  rates  (960/s  to  46,400/s). 


This  GLO  process  was  repeated  a  number  of  times,  with  material  constants  chosen  by  GLO  but 
guided  by  constraints  set  by  the  user  (maximum,  minimum,  starting  value,  etc.),  until  a  minimum 
FOM  is  reached.  After  this  automated  GLO  process  was  completed,  a  manual  search  was  done 
in  the  vicinity  of  the  minimum  found  by  GLO  to  further  refine  the  solution.  Improvement  was 
small  but  useful. 
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One  final  note  is  worth  mentioning.  Examination  of  equation  9  shows  that  for  initial  yield 
(plastic  strain  =  0),  the  third  term  is  always  zero.  Thus,  early  attempts  divided  the  optimization 
into  fitting  the  remaining  constants  (Co,  C\,  (3o,  and  sp )  to  just  the  initial  yield  information  from 

the  five  data  sets,  then  holding  them  fixed  while  fitting  the  strain  hardening  constants.  This 
approach  yielded  much  higher  (worse)  figures  of  merit  than  a  global  approach,  wherein  all 
parameters  were  allowed  to  vary  at  once.  The  reason  for  this  is  that  the  optimization  process  was 
constrained  by  the  (albeit  accurate)  fit  to  the  initial  yield  data  points,  and  GLO  was  not  allowed 
to  vary  those  few  points  in  order  to  improve  the  overall  fit  at  all  strains.  Since  the  high  strain 
performance  of  the  model  is  important  in  ballistic  events,  a  global  approach  was  used. 


4.  The  Fit 


The  optimization  process  yielded  the  set  of  material  constants  in  table  1  for  equation  9,  for  the 
low  cost  Ti-6A1-4V  material. 


Table  1.  Parameters  for  equation  9  for  low  cost  Ti-6A1-4V. 


C0,  MPa 

1217 

C2,  MPa 

3955 

Cb  MPa 

139 

sr 

0.1877 

£ p  ,/S 

1.597e4 

£(x  ^ 

8.249e5 

Po,/K 

1.591e-2 

a0,  /K 

7.549e-3 

The  FOM  for  this  fit  was  0.0107  or  an  average  error  of  only  1.07%.  The  fit  obtained  by  GLO 
when  an  exponent  of  1/2  was  used  (as  in  equation  3)  resulted  in  a  FOM  of  0.015 1  or  about  40% 
worse  than  that  obtained  with  an  exponent  of  1 .  Importantly,  a  value  of  1/2  led  to  lower  values 
of  ea  and  sp ,  which  is  detrimental  to  the  ability  of  the  model  to  predict  localization,  as 

previously  discussed  (and  a  quantitative  example  of  this  is  given  in  a  later  section).  Plotting  the 
strain  function  in  equation  9  (i.e.,  exponent  =  1)  shows  that  saturation  of  the  function  occurs  at  a 
strain  of  about  0.70  for  er  =  0.20  and  at  a  strain  of  about  2.00  for  er  =  0.50. 

The  fit  compared  to  the  discretized  data  is  shown  in  figure  1 .  The  fit  to  the  0. 1/s  data  (which  was 
omitted  from  the  FOM  calculations  as  well  as  the  calibration  process)  is  shown  for  information. 
Even  though  0.1/s  was  omitted  from  the  fitting  process,  the  model  prediction  still  passes  through 
some  of  the  0.1/s  data  points.  This  indicates  that  the  model’s  predictions  of  strain  rate  hardening 
and  thermal  coupling  (since  adiabatic  heating,  important  at  the  high  rates,  is  negligible  at  the  low 
rate)  are  well  represented  over  the  range  of  the  data. 
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The  strain  rate  effect  is  shown  in  figure  2;  the  bi-linear  aspect  of  the  data  is  apparent.  The  9585/s 
data  set  is  near  the  critical  point  of  the  plot.  Since  equation  9  is  a  continuous  function,  it  is  not 
able  to  closely  represent  the  sharp  cusp7  in  the  stress  in  the  vicinity  of  this  strain  rate.  This 
difficulty  is  apparent  in  figure  1,  where  the  fit  to  the  9585/s  curve  shows  the  most  error  of  all  the 
high  rates.  As  if  to  emphasize  the  critical  nature  of  this  strain  rate,  even  the  data  show  a  high 
degree  of  variability.  The  maximum  and  minimum  individual  tests  are  shown  as  dashed  lines  in 
figure  1 ,  and  the  minimum  actually  overlaps  some  lower  rate  data.  The  model  prediction  falls 
within  these  limits  of  the  experimental  variation.  The  prediction  of  initial  yield  for  the  0.1 /s  case 
is  in  error  by  about  15%  in  figure  2,  but  this  is  not  considered  a  liability  for  ballistic  applications; 
sufficient  accuracy  is  attained  by  10  /s,  and  the  important  processes  in  ballistic  problems 
typically  occur  at  strain  rates  greater  than  that. 
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Figure  1.  Modified  Zerilli- Armstrong  model  predictions  compared  to  experimental  data. 


The  quality  of  the  fit  to  the  initial  yield  strength  can  be  examined  in  another  way  if  predicted  initial 
yield  stress  is  plotted  against  the  exponent  of  the  second  term  in  equation  9  as  the  independent 
variable.  This  exponent  provides  the  functional  dependence  of  initial  yield  on  temperature  and 
strain  rate  (recall  that  the  third  term  does  not  play  a  role  in  initial  yield  strength),  but  Casern’s  data 
are  at  room  temperature  (293  K,  i.e.,  no  temperature  preconditioning  and  no  adiabatic  heating  for 
initial  yield).  Figure  3  shows  the  comparison;  high  rate  data  are  at  the  left.  The  model  captures  the 
data  nicely,  with  the  exception  of  the  0.1/s  experiment,  seen  at  the  right  in  the  plot. 


7The  cusp  is  probably  continuous  in  the  vicinity  of  9585/s,  although  the  slope  is  changing  rapidly.  The  curve  in 
figure  2  is  drawn  as  discontinuous  for  simplicity. 
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Figure  2.  Dependence  of  initial  yield  on  the  strain  rate. 
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Figure  3.  Evaluation  of  the  prediction  of  the  initial  yield  stress  based  on  the  form  of  the  model. 


For  comparison,  the  same  GLO  method  was  used  to  fit  the  Johnson-Cook  model  to  Casern’s 
data.  Table  2  gives  the  resulting  material  constants  for  the  function  in  equation  1.  The  figure  of 
merit  for  this  fit  was  0.023  or  an  average  error  of  2.3%.  Figure  4  shows  the  performance  of  the 
model  compared  to  the  experimental  data.  The  Johnson-Cook  model  allows  only  a  single  value 
for  the  material  constant  C  in  equation  1 .  Since  C  is  proportional  to  the  slope  in  figure  2,  the 
single  “average”  slope  of  the  model  allows  reasonable  predictions  for  three  strain  rate  curves  but 
misses  two  of  the  high  strain  rate  curves.  Modifying  the  model  to  use  a  bi-linear  slope  as  in 
figure  2  improved  the  fit  to  the  Hopkinson  Bar  data  but  provided  no  improvement  in  predictions 
of  a  Taylor  anvil  test.  As  with  the  modified  Zerilli- Armstrong  model,  the  Johnson-Cook  model 
also  has  trouble  with  the  critical  9585/s  rate.  Examination  of  figure  4  also  shows  that  the 
Johnson-Cook  model  does  not  capture  stress  saturation  at  high  strain  as  well  as  the  modified 
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Zerilli- Armstrong  model.  Interestingly,  n  =  0.812,  close  to  the  exponent  of  1  for  the  strain 
function  in  equation  9. 


Table  2.  Parameters  for  equation  1 
for  low  cost  Ti-6A1-4V. 


A,  MPa 

840 

B,  MPa 

550 

C 

0.0664 

n 

0.812 
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1.769 
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Figure  4.  Johnson-Cook  predictions  compared  to  experimental  data. 


5.  Testing  the  Fit  in  Simulations 


Equation  9  and  parameters  from  table  1  were  tested  in  a  series  of  simulations  with  axisymmetric 
EPIC.  Simulations  of  four  sets  of  experiments  were  conducted:  Taylor  anvil  experiments,  V50 
experiments  using  a  steel  projectile,  penetration  of  a  Ti-6A1-4V  penetrator  into  semi-infinite  Ti- 
6A1-4V,  and  penetration  of  a  tungsten  penetrator  into  semi-infinite  Ti-6A1-4V.  All  simulations 
used  the  Mie  Gruneisen  equation  of  state,  with  EPIC  library  constants,  for  all  materials.  All 
simulations  (except  some  simulations  of  the  Taylor  anvil  experiments;  see  the  related  footnote  in 
the  next  section)  used  the  Johnson-Cook  fracture  model,  with  EPIC  library  constants,  for  all 
materials.  All  simulations  used  EPIC’s  General  Particle  Algorithm,  except  for  the  Taylor  anvil 
simulations. 
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5.1  Simulation  of  Taylor  Anvil  Experiments 

In  addition  to  the  Hopkinson  Bar  experiments,  Casern  conducted  Taylor  anvil  experiments  using 
specimens  fabricated  from  the  low  cost  Ti-6A1-4V  material  (14).  A  summary  of  the  results  of 
several  experiments  at  205  m/s  and  axisymmetric  EPIC  simulations  (10  crossed  elements  across 
the  radius  of  the  specimen)  with  the  modified  Zerilli- Armstrong  model  and  the  Johnson-Cook 
model  is  given  in  table  3.  Overall,  both  models  accurately  described  the  final  deformed  shape  of 
the  specimen.  Also  included  in  the  table  is  the  result  of  a  simulation  using  the  EPIC  library 
parameters  for  Ti-6A1-4V  in  the  Johnson-Cook  model.  Those  parameters  were  not  intended  for 
the  low  cost  material  and  are  included  to  satisfy  curiosity.  Actually,  the  result  is  fairly  good, 
although  somewhat  worse  than  the  models  that  were  calibrated  specifically  for  this  material. 
Finally,  for  completeness,  the  prediction  of  the  Zerilli- Armstrong  bcc  model  from  reference  7  is 
included.  This  bcc  model  was  calibrated  with  a  different  set  of  Hopkinson  Bar  data8,  one  that 
included  strain  rates  from  quasi-static  to  only  3000/s.  Overall,  it  was  slightly  less  accurate  than 
the  present  modified  Zerilli-Armstrong  model  for  replicating  this  experiment. 

Table  3.  Summary  of  the  Taylor  anvil  experiments  and  simulations. 


Experiment 

Modified 

Zerilli-Armstrong 

Calibrated 

Johnson-Cook 

EPIC 

Johnson-Cook 

bcc 

Zerilli-Armstrong 

mm 

mm 

percent  error 

mm 

percent  error 

mm  |percent  error 

mm 

percent  error 

Final  Length 

16.93 

16.95 

0.1 

16.97 

0.2 

16.31  1  3.8 

16.64 

1  L7 

Maximum  Diameter 

3.62 

3.54 

2.2 

3.57 

1.4 

3.94  |  8.8 

3.56 

1  1.7 

More  details  of  the  EPIC  simulation  of  the  Taylor  anvil  test  using  the  modified  Zerilli-Armstrong 
model  are  given  in  figure  5.  On  the  left  is  a  view  of  the  strain  field  at  14  ps,  just  after  the  specimen 
separates  (bounces)  from  the  anvil.  In  the  center  is  the  profile  of  the  specimen  at  14  ps.  Experi¬ 
mental  data  points  are  plotted  on  the  figure,  showing  quite  good  agreement  between  simulation  and 
experiment.  The  simulation  shows  a  more  pronounced  effect  of  the  friction  at  the  interface  between 
the  specimen  and  the  anvil.  On  the  right  in  figure  5  is  a  plot  of  strain  rate  at  1  ps  when  strain  rates 
are  generally  the  highest.  It  shows  that  most  of  the  impact  end  of  the  specimen  is  strained  at  rates 
of  104/s  to  105/s.  However,  these  high  rates  are  short  lived,  dropping  to  a  maximum  of  2xl04/s  by 
5  ps  and  continuing  to  drop  rapidly. 

Some  but  not  all  of  the  experiments  at  205  m/s  (nominally)  showed  shear  localization  in  the 
specimen;  figure  6  shows  a  typical  example.  Apparently,  this  impact  velocity  is  a  marginal 
condition  for  causing  this  type  of  failure.  (The  experimental  results  listed  in  table  3  are  average 
values  for  two  specimens  that  did  not  localize.) 


O 

The  bcc  model  was  not  recalibrated  with  the  present  Hopkinson  Bar  data  because  equation  9  was  developed  for 
hep  Ti-6A1-4V  and  supersedes  the  bcc  model  for  this  material.  Alternatively,  the  Johnson-Cook  model  makes  no 
such  distinction  (i.e.,  it  should  be  applicable  to  Ti-6A1-4V),  so  it  was  calibrated  here  for  comparison  purposes. 
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Figure  5.  Details  of  the  EPIC  simulation  of  the  Taylor  anvil  test  with  the  use  of  the  modified 
Zerilli-Armstrong  model. 


Figure  6.  Photograph  of  shear  localized  Taylor  anvil  specimen  (impact  velocity  of  205  m/s). 


The  modified  Zerilli-Armstrong  model  did  not  predict  localization  for  this  impact  velocity  with 
the  parameter  set  of  table  2  (see  figure  5)9.  To  investigate  what  would  be  required  for  the  model 
to  predict  localization  in  the  simulation,  two  excursions  were  undertaken.  In  the  first,  the  impact 
velocity  was  increased  in  several  sequential  steps.  Figure  7  shows  that  localization  appears 


9The  prediction  of  localization  in  the  EPIC  simulations  of  figures  7  and  8  required  FAIL=1,  whereas  the  simulation 
of  figure  5  used  FAIL=0  (it  was  the  only  simulation  reported  here  that  used  FAIL=0).  FAIL=0  will  not  allow  fracture 
of  the  material  when  damage  exceeds  1  in  the  Johnson-Cook  fracture  model,  and  FAIL=1  allows  degrading  of  the 
material  strength  when  damage  exceeds  1 . 
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abruptly  at  221  m/s,  along  the  plane  of  failure10  seen  in  figure  6.  This  localization  begins  when 
the  predominant  strain  rate  drops  from  a  maximum  of  about  6xl04/s  to  under  2xl04/s,  roughly 
equivalent  to  the  value  of  sp ,  the  lower  of  the  two  strain  rate  constants  (table  1).  This  occurs 

roughly  9  ps  after  impact. 


220m/s 


221m/s 


CONTOURS 

0.00E+00 

■  1.22E-01 

■  2.44E-01 
3.66E-01 

■  4.88E-01 

■  5.80E-01 


Figure  7.  Predicted  strain  contours  showing  incipient  strain  localization  in  the  Taylor  anvil  test  (t  =  1 1  ps  is 
shown;  separation  from  the  anvil  occurs  at  abour  14  ps). 


At  slightly  higher  velocity,  the  localization  seen  in  figure  7  extends  farther  from  the  impact  face, 
but  the  220-m/s  case  did  not  show  localization  at  any  time  (the  simulation  was  run  until  separation 
of  the  specimen  from  the  anvil,  15  ps  after  impact).  The  model  prediction  is  consistent  with  two 
experimental  observations.  First  is  the  location  and  orientation  of  the  localization  as  seen  in 
figures  6  and  7.  Second  is  the  incipient  nature  of  the  localization  since  the  model  is  reasonably 
predictive  of  its  onset  (predicting  the  onset  velocity  within  about  8%  of  the  experiment).  The 
important  achievement  is  the  actual  prediction  of  localization  by  the  strength  model. 

Variations  of  Po  were  tested  in  the  second  excursion  to  investigate  prediction  of  localization.  Each 
variation  of  po  was  balanced  by  adjustments  in  sp  in  order  to  keep  Pi  unchanged  (in  equation  6).  If 

Po  is  increased  by  just  20%,  (po  =  1.909e-2,  sp  =  1.105e5/s),  localization  is  observed  in  the  simu¬ 
lation  (see  figure  8).  The  resulting  sp  is  on  the  order  of  the  highest  strain  rates  seen  in  the  Taylor 

anvil  simulation  of  figure  5.  This  small  change  in  constants  is  also  consistent  with  a  marginal 
localization  condition  seen  in  the  experiments.  If  Po  is  increased  by  50%,  (Po  =  2.387e-2, 

Sp  =  2.022e6/s),  localization  is  more  organized  and  occurs  at  an  earlier  time,  as  seen  in  figure  8.  In 

this  case,  sp  is  greater  than  the  highest  strain  rates  seen  in  figure  5,  and  this  is  most  likely  the  cause 

of  the  improved  localization.  The  deformation  of  the  specimen  seemed  to  be  ending  at  about  9.5  ps 
(see  the  inset  graph  in  figure  8).  When  the  localization  formed  at  1 1  ps,  deformation  continued  in 
earnest.  The  simulations  of  figure  8  are  the  only  uses  of  these  sets  of  constants  reported  here. 

Subsequently,  the  optimization  of  the  model  constants  was  re-run  with  this  latter  value  of  Po  and 
Sp  fixed  (FOM  was  1.67%).  A  simulation  using  this  set  of  constants  did  not  show  localization, 
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The  term  “plane”  is  used  loosely  in  the  context  of  the  simulation  since  the  simulation  is  an  axisymmetric 


calculation. 
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most  likely  because  of  the  reduced  value  of  sa  =  5.13e4  that  resulted  from  the  GLO  optimiza-tion. 

This  low  value  (relative  to  8.249e5,  the  value  from  table  1  used  for  figure  8)  overpowers  the 
increased  value  of  ep  and  allows  rate  hardening  to  overtake  thermal  softening  too  soon. 

Prediction  of  localization  thus  requires  both  a  and  p  to  be  positive  at  the  highest  strain  rates  seen 
in  the  problem,  or  at  least  the  combined  effect  must  be  positive  (e.g.,  one  can  be  slightly  negative). 


Figure  8.  Localization  in  the  Taylor  anvil  simulation  attained  when  the  value  of  p0  and  £  „ 
was  increased. 

5.2  Simulation  of  a  V50  Experiment 

Burkins  et  al.  (11)  studied  the  effect  of  annealing  temperature  on  the  ballistic  limit  velocity  of 
28.5-mm-thick  Ti-6A1-4V  plates  attacked  by  a  20-mm  steel  fragment  simulating  projectile 
(FSP).  They  found  the  V50  of  this  configuration  to  be  between  1092  m/s  and  1 141  m/s  for  the 
various  annealing  conditions  studied.  In  those  shots  that  resulted  in  perforation  of  the  target, 
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failure  of  the  target  was  observed  to  occur  by  plugging  (i.e.,  shearing  out  of  a  cylindrical  plug  of 
target  material  ahead  of  the  penetrator).  In  those  shots  that  resulted  in  incomplete  penetration, 
shear  bands  were  often  observed  parallel  and  perpendicular  to  the  direction  of  fire.  Figure  9 
shows  a  typical  example  (from  reference  11;  reprinted  by  permission  of  the  principal  author). 


Figure  9.  Sectioned  view  of  the  impact  crater  and  remaining  plate  from  a  V50  experiment  (taken  from 
reference  (77);  reprinted  by  permission).  (Striking  velocity  was  1 106  m/s,  near  the  limit 
velocity.) 


The  test  configuration  was  simulated  in  axisymmetric  EPIC  with  the  use  of  a  20-mm-long  by 
20-mm-diameter  steel  cylinder  (10  crossed  elements  across  the  FSP  radius)  and  the  model  of 
equation  9  and  table  1.  Results  of  two  significant  simulations  are  shown  in  figure  10.  At 
1020  m/s,  the  FSP  was  unable  to  perforate  the  target  (notice  that  it  has  begun  to  rebound  out  of 
the  target  at  120  ps),  although  some  spalling  occurs  on  the  rear  surface.  At  1040  m/s,  spalling 
also  occurs  at  the  early  time,  but  perforation  is  clearly  occurring  at  the  late  time.  Thus,  a 
predicted  V5o  of  approximately  1030  m/s  is  indicated,  which  is  within  6%  to  11%  of  the 
experimental  values. 

The  most  striking  feature  in  figure  10  is  the  appearance  of  strain  localization  at  the  1040-m/s 
impact  velocity.  Thus,  the  threshhold  perforation  at  1040  m/s  is  accomplished  by  strain  localiza¬ 
tion  in  the  target  and  plugging  of  the  rear  portion  of  the  target,  consistent  with  what  was  observed 
experimentally  by  Burkins  et  al.  (11,  15).  The  simulation  even  shows  a  small  amount  of  the 
perpendicular  shear  banding,  which  indicates  that  the  origin  of  this  failure  seen  in  the  experiments 
(figure  9)  is  more  than  just  metallurgical  (e.g.,  from  the  rolling  of  the  plates  during  processing). 
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Strain  rate  contours  in  figure  10  show  that  the  highest  rates  in  the  1040-m/s  case  are  about  105/s,  a 
strain  rate  slightly  greater  than  sp  and  slightly  less  than  sa  (table  1). 
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Figure  10.  Localization  and  plugging  in  simulations  of  V50  experiments  with  the  use  of  equation  9  to  model 
the  target.  (Experimental  V50  for  this  target  was  approximately  1100  m/s.) 


When  the  impact  velocity  was  increased  to  1300  m/s,  a  second  cylindrical  localization  band 
appeared  at  a  slightly  larger  diameter  than  the  primary  band  (that  of  figure  10).  The  slight  per¬ 
pendicular  shear  banding  observable  in  figure  10  for  the  1040-m/s  case  became  slightly  more 
pronounced.  At  1400  m/s,  the  localization  and  plugging  continued  to  mark  the  failure  process, 

5  5  and  the  small  perpendicular  shear  band  became  slightly  more  pronounced.  The  plugging  and 
localization  continued  until  at  least  1500  m/s,  but  at  a  velocity  of  2000  m/s,  no  plugging  or 
localization  was  observed — just  massive  deformation  and  failure  as  the  FSP  pushed  violently 
through  the  target. 

For  comparison,  the  Burkins’  et  al.  V5o  experiment  was  modeled  with  the  Zerilli-Armstrong  bcc 
model  from  reference  7.  The  V50  for  the  bcc  model  was  about  1400  m/s,  a  significant  error. 
Even  though  that  model  did  not  include  any  explicit  measures  to  model  localization  as  does 
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equation  9,  the  bee  model  did  show  the  beginning  of  localization.  Figure  1 1  shows  a  good 
example.  However,  the  localization  seen  in  figure  1 1  did  not  grow  at  times  beyond  40  ps,  and  at 
higher  velocity,  the  localization  was  swamped  by  widespread  failure  of  the  material  ahead  of  the 
penetrator.  Thus,  the  localization  did  not  develop  into  plugging  as  seen  in  the  experiment  (figure 
9)  and  in  the  simulation  (figure  10).  One  reasonable  interpretation  of  this  outcome  is  that  the 
coupling  of  strain  rate  hardening  and  thermal  softening  is  well  represented  in  both  the  bee  model 
and  the  present  model  (the  coupling  is  in  the  exponential  terms,  which  are  similar  in  the  two 
models),  but  the  Taylor  strain  hardening  of  the  bee  model  inhibited  growth  of  the  localization, 
whereas  the  strain  function  of  the  present  model  favors  that  growth. 
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Figure  11.  Nascent  localization  with  the  Zerilli- Armstrong  bcc  model  in  simulations  of  V50  experiments 
(1200-m/s  impact  velocity  shown). 


5.3  Simulation  of  Ti-6A1-4V  on  Ti-6A1-4V  Penetration  Experiments 

Meyer  and  Kleponis  (7)  performed  penetration  experiments  using  semi-infinite  targets  fabricated 
from  the  Ti-6A1-4V  material,  impacted  by  penetrators  made  of  the  same  material.  A  summary  of 
a  portion  of  their  results  is  shown  in  figure  12.  Their  experimental  results  are  shown  as  the  solid 
symbols,  and  the  line  is  a  fit  to  the  data.  The  open  triangles  show  the  results  of  3-D  CTH  simula¬ 
tions  using  the  Zerilli- Armstrong  bcc  model,  included  here  for  comparison.  The  open  squares 
show  the  results  for  axisymmetric  EPIC  (five  crossed  elements  across  the  penetrator  radius)  with 
the  current  model.  Predicted  depth  of  penetration  (DOP)  is  reasonably  accurate  from  1200  m/s  to 
1600  m/s,  but  the  error  is  substantial  at  2000  m/s.  Maximum  strain  rates  at  this  impact  velocity  are 
greater  than  106  in  the  target,  higher  than  the  data  used  to  calibrate  the  model,  and  105  in  the 
penetrator,  slightly  higher  than  the  data.  Maximum  rates  were  about  one-third  as  high  at  an  impact 
velocity  of  1600  m/s.  Efforts  are  under  way  to  install  this  model  into  CTH  in  order  to  conduct  3-D 
simulations  for  direct  comparison  to  the  bcc  model  predictions. 
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Figure  12.  Simulation  results  compared  to  experiment  and  3-D  CTH  with  an  earlier  ZA  model. 


5.4  Simulation  of  the  Penetration  of  the  190W  Round  Into  Ti-6A1-4V 

Gooch  and  Burkins  ( 16)  studied  the  performance  of  Ti-6A1-4V  attacked  by  various  penetrators, 
including  the  190W,  which  is  an  l/d=20,  93%  tungsten  rod,  190.5  mm  long.  Impact  velocities 
for  the  190W  were  between  950  m/s  and  1700  m/s.  Gooch  and  Burkins  recorded  the  DOPs  and 
fit  those  data  with  an  analytical  function,  plotted  as  the  solid  line  in  figure  13.  Plotted  as  open 
symbols  in  figure  13  are  results  of  axisymmetric  EPIC  simulations  (five  crossed  elements  across 
the  penetrator  radius)  with  three  different  strength  models  for  the  Ti-6A1-4V.  The  190W  was 
modeled  with  the  Johnson-Cook  strength  model  using  constants  from  reference  4  (except  the 
yield  strength  here  was  1.17  GPa),  the  Mie-Gruneisen  equation  of  state  (EOS),  and  the  Johnson- 
Cook  fracture  model.  EPIC  library  constants  were  used  for  the  EOS  and  fracture  model. 


Figure  13.  Axisymmetric  EPIC  simulation  results  for  the  current  model  compared  to  experiment  and 
other  models. 


Predictions  of  depth  of  penetration  into  the  Ti-6A1-4V  were  inaccurate  when  the  Ti-6A1-4V  was 
modeled  with  the  Johnson-Cook  strength  model  using  constants  from  the  EPIC  library  (open 
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circles  in  figure  13).  This  shows  the  error  of  using  that  constant  set  to  model  the  low  cost  Ti- 
6A1-4V  material.  The  300-mm-thick  target  was  perforated  at  1700  m/s  with  that  model.  When 
the  set  of  constants  in  table  2  is  used  with  the  Johnson-Cook  strength  model,  results  are  much 
improved  (open  triangles  in  figure  13),  although  inaccuracy  increases  as  impact  velocity  is 
increased.  Results  are  much  better  over  the  range  of  1 100  m/s  to  2000  m/s  when  the  modified 
Zerilli- Armstrong  model  is  used  with  constants  from  table  1 . 


6.  Conclusion 


The  Zerilli-Armstrong  hep  model  was  modified  and  fit  to  high  rate  Hopkinson  Bar  data.  The  fit 
to  the  data  was  very  good.  The  model  was  installed  into  EPIC  and  performed  well  in  a  varied 
series  of  axisymmetric  simulations.  It  did  not  perform  well  for  the  case  of  a  Ti-6A1-4V  penetrator 
impacting  a  block  of  Ti-6Al-4  V  at  2000  m/s. 


7.  Recommendation 


Because  CTH  is  used  for  large  scale  penetration  simulations  at  ARL,  the  model  should  be  installed 
into  CTH  and  tested  in  the  same  series  of  simulations  with  3-D  CTH.  Numerical  effects  of  the 
axisymmetric  calculations  (e.g.,  anomalies  on  the  axis  of  symmetry)  in  EPIC  may  have  affected 
the  high-velocity  EPIC  penetration  simulations,  and  this  provides  additional  motivation  to  examine 
the  model’s  performance  in  3-D  CTH. 
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